Análisis Causal de la Mortalidad de Renacuajos de Rana Carrizo

Métodos Analíticos

Autores/as

Blanca Estela García Manjarrez - 118886

Yuneri Pérez Arellano – 199813

Fecha de publicación

21 de mayo de 2025

Introducción

En este proyecto empleamos los datos experimentales de Vonesh & Bolker (2005)1, quienes en su investigación examinaron las consecuencias de la plasticidad de eclosión inducida por depredadores desde la etapa larval hasta la metamorfosis en la rana de caña de África Oriental, Hyperolius spinigularis2 realizando un experimento en el que manipulaban el tamaño y la densidad larvaria inicial (imitando los efectos de los depredadores de los huevos). Esperaban que las crías inducidas por depredadores (porque están menos desarrolladas y son más pequeñas) experimentaran mayores tasas de depredación per cápita y un período larvario más largo y, por lo tanto, exhibirían una menor supervivencia a la metamorfosis en presencia de depredadores acuáticos que las larvas más grandes, más desarrolladas y eclosionadas más tarde. Sin embargo, los resultados mostraron que las larvas inducidas por depredadores no solo sobrevivieron a la metamorfosis, sino que también tuvieron tasas de crecimiento más rápidas y alcanzaron tamaños más grandes en la metamorfosis. Esto los motivó a desarrollar un modelo parametrizado a partir de experimentos adicionales para explorar si una combinación de mecanismos, crecimiento compensatorio y depredación específica por densidad y tamaño, podría dar lugar a este patrón. Es por eso que con esta introducción, buscamos replicar y entender su trabajo, utilizando un enfoque bayesiano jerárquico para modelar la mortalidad larval y las respuestas compensatorias posteriores. En este sentido, el modelo jerárquico nos permitirá capturar la heterogeneidad entre los tanques de renacuajos y compartir información entre ellos, lo que resulta en estimaciones más robustas y precisas.

Datos

Los datos provienen de la librería de rethinking de R. Constan de 48 observaciones, que representan los tanques de renacuajos clasificados en pequeños, medianos y grandes, dependiendo de la densidad de renacuajos en cada uno. Además, de información sobre la supervivencia (variable binaria) y de la tasa de supervivencia en cada tanque.

Variable Descripción
density Densidad inicial de renacuajos
pred Factor indicador de presencia de depredadores
size Tamaño de los renacuajos
surv Número de renacuajos que sobrevivieron
propsurv Proporción de supervivencia (surv/density)

En este experimento observamos mucha variación en los datos, y no toda se debe al tratamiento experimental (como la presencia de depredadores). Una gran parte de esa variación proviene de factores no medidos, propios del entorno donde viven los renacuajos. Podemos imaginar cada fila del dataset, como un “tanque”, es decir, un pequeño ambiente experimental que contiene varios renacuajos. Aunque algunos tanques tengan la misma densidad o condiciones aparentes, hay muchas cosas que no estamos midiendo (como temperatura, luz, microalgas, etc.) que también influyen en la tasa de supervivencia. Esto hace que los tanques funcionen como lo que llamamos un conglomerado o cluster. Dentro de cada tanque observamos múltiples renacuajos, por lo que los datos tienen una estructura agrupada. En otras palabras, tenemos medidas repetidas dentro de grupos que son diferentes entre sí.

Para nuestro análisis, nos centraremos en surv como variable de respuesta (binomial) frente a density como total de ensayos.

Es importante mencionar, que si usamos el mismo valor base (intercepto) para todos los tanques (pooling), estamos ignorando diferencias importantes entre ellos. Esto puede hacer que no detectemos correctamente el efecto de otras variables como la densidad o el predador. Si por el contrario, usamos un intercepto distinto para cada tanque (no pooling), pero sin compartir información entre ellos, podríamos caer en lo que se llama una “amnesia estadística”: tratando a cada tanque como si no se tuviera nada que aprender de los demás. Pero eso no tiene sentido, porque aunque cada tanque es diferente, los datos de uno pueden ayudarnos a entender mejor a los demás.

Por ello, empleamos también un modelo bayesiano jerárquico o multinivel o como lo mencionaremos en este proyecto: modelo con interceptos variables (partial pooling), de este modo, cada tanque tiene su propio parámetro de línea base, y al mismo tiempo estimamos la dispersión entre tanques mediante un prior adaptativo, que aprende de los datos. Con esto buscamos lograr un equilibrio, entre asumir que todo es igual (subajuste) y asumir que todo es completamente distinto (sobreajuste).

Los objetivos de este proyecto son:

  1. Reproducir y comprender el ejemplo de Statistical Rethinking aplicado a los datos de Reedfrogs.
  2. Modelar la mortalidad larval y las respuestas compensatorias posteriores explorando distintos niveles de agrupamiento pooling, no pooling y partial pooling.
  3. Explorar el trade-off underfitting/overfitting comparando los tres niveles de pooling (pooling, no pooling y partial pooling) mediante criterios de información y gráficos de predicción.
  4. Desplegar un análisis causal formal con un DAG que recoja nuestros supuestos de identificación.

Con este enfoque se busca profundizar en los costes y beneficios de la eclosión temprana inducida por depredadores, y demostrar cómo la regularización adaptativa de los modelos jerárquicos, permite inferir efectos individuales de forma más robusta, en presencia de datos jerarquizados y dispersos.

DAG

A continuación, se presenta la Gráfica Dirigida Acíclica (DAG) que ilustra las relaciones causales entre las variables de interés. Este DAG se basa en la premisa de que la densidad de renacuajos, el tamaño y la presencia de depredadores influyen en la supervivencia, y que la jerarquía de los tanques también afecta a estas relaciones.

DAG causal: densidad, tamaño, depredadores y jerarquía de tanque

Con

  • \(\textrm{T}=\textrm{Tanque}\)

  • \(\textrm{D}=\textrm{Densidad inicial}\)

  • \(\textrm{G}=\textrm{Tamaño}\)

  • \(\textrm{P}=\textrm{Depredadores}\)

  • \(\textrm{S}=\textrm{Supervivencia}\)

Aunque podríamos imaginar otras dependencias entre las variables, hay que tener presente que estos datos provienen de un experimento controlado, es decir, como se expone en la introducción, el experimento es manipulado y cerrado bajo las condiciones que Vonesh & Bolker establecieron. En un escenario natural, sería razonable investigar vínculos, como el efecto del tamaño de los renacuajos en la densidad poblacional, la influencia de los depredadores sobre esa densidad, o el papel de variables no registradas —por ejemplo, la disponibilidad de alimento u otros recursos— tanto en el tamaño como en la densidad, e incluso factores genéticos que modulen el desarrollo de los renacuajos. Sería muy valioso, repetir estas estimaciones en cuerpos de agua naturales, en lugar de tanques de laboratorio. Por ahora, este experimento, nos permite centrarnos en la tasa de supervivencia, bajo condiciones estrictamente controladas, estableciendo una base sólida para futuros estudios en la naturaleza.

Modelos

1. Modelo Total Pooling

En este primer modelo totalmente agrupado, se asume que todos los tanques tienen la misma probabilidad de supervivencia. No hay diferencias entre tanques, salvo la variación por la densidad inicial \(D_i\).

\[ S_i \sim \textrm{Binomial}(D_i,p_i) \]

\[ \textrm{logit}(p_i) = \alpha \]

\[ \alpha = \textrm{Normal}(0, 1.5) \] Este modelo ignora la heterogeneidad entre tanques (total pooling) y servirá como línea base para comparar con el modelo jerárquico.

data {
  int<lower=0> T;         // Num de tanques
  int<lower=0> S[T];         // Num de renacuajos que sobrevivieron
  int<lower=0> D[T];         // Densidad inicial
}

parameters {
  real alpha;                // Un alpha para todos los tanques
}

model {
  // Prior de aplha
  alpha ~ normal(0, 1.5);
  
  for (t in 1:T) {
    S[t] ~ binomial(D[t], inv_logit(alpha));
  }
}

generated quantities {
  int S_rep[T]; 

  // Predicciones basadas en probabilidad comun
  for (t in 1:T) {
    S_rep[t] = binomial_rng(D[t], inv_logit(alpha));
  }
}

En el enfoque de total pooling, asignamos una única \(\alpha\) a todos los tanques, de modo que el modelo asume idéntica probabilidad de supervivencia para los renacuajos en cada uno de ellos. Es como si tratáramos todos los tanques como réplicas exactas, sin captar ninguna heterogeneidad más allá de la variación provocada por la densidad inicial.

A continuación, mostramos los resultados obtenidos para la estimación de \(\alpha\) en nuestro modelo inicial.

parameter mean sd 2.5% 50% 97.5% n_eff Rhat
alpha 0.84 0.06 0.72 0.84 0.97 1,255.30 1.00
lp__ −685.65 0.66 −687.58 −685.40 −685.17 2,318.65 1.00

Podemos cotejar nuestras predicciones con los datos originales:

Predicciones del modelo 1 (total pooling)

En este modelo totalmente agrupado (total pooling), las predicciones (puntos rojos) se limitan a la media global de supervivencia (línea punteada azul) y no siguen la dispersión real de los tanques (puntos grises), lo que indica un claro subajuste.

2. Modelo No-Pooling

En este modelo no agrupado (no pooling) asignamos un intercepto \(\alpha_i\) distinto a cada tanque, pero no compartimos información entre ellos. Esto significa que cada tanque tiene su propio parámetro de línea base, y no hay aprendizaje entre ellos.

\[ S_i \sim \textrm{Binomial}(D_i,p_i) \]

\[ \textrm{logit}(p_i) = \alpha_i \]

\[ \alpha_i = \textrm{Normal}(0, 1.5)\]

Lo cual se traduce en el siguiente código de Stan:

data {
  int<lower=0> T;             // Número de tanques
  int<lower=0> S[T];          // Número de sobrevivientes
  int<lower=0> D[T];          // Densidad inicial
}

parameters {
  real alpha[T];              // alpha para cada tanque
}

model {
  // Priors para cada alpha_i
  for (i in 1:T) {
    alpha[i] ~ normal(0, 1.5);
  }
  
  // Modelo para cada tanque
  for (t in 1:T) {
    S[t] ~ binomial(D[t], inv_logit(alpha[t]));
  }
}

generated quantities {
  int S_rep[T];

  // Generar predicciones 
  for (t in 1:T) {
    S_rep[t] = binomial_rng(D[t], inv_logit(alpha[t]));
  }
}

Podemos ver los resultados de el modelo No-Pooling para la estimación de \(\alpha\).

parameter mean sd 2.5% 50% 97.5% n_eff Rhat
alpha[1] 1.71 0.77 0.32 1.67 3.38 4,727.77 1.00
alpha[2] 2.40 0.89 0.85 2.35 4.29 4,651.05 1.00
alpha[3] 0.76 0.66 −0.46 0.74 2.12 5,454.85 1.00
alpha[4] 2.41 0.89 0.80 2.37 4.33 5,199.02 1.00
alpha[5] 1.71 0.76 0.29 1.67 3.30 4,768.40 1.00
alpha[6] 1.71 0.76 0.35 1.66 3.37 4,306.38 1.00
alpha[7] 2.42 0.90 0.83 2.36 4.41 4,062.44 1.00
alpha[8] 1.71 0.74 0.38 1.66 3.26 5,347.54 1.00
alpha[9] −0.38 0.61 −1.59 −0.37 0.81 4,527.06 1.00
alpha[10] 1.72 0.77 0.35 1.66 3.36 4,798.29 1.00
alpha[11] 0.75 0.64 −0.44 0.73 2.04 4,876.00 1.00
alpha[12] 0.36 0.62 −0.81 0.36 1.59 4,630.22 1.00
alpha[13] 0.75 0.65 −0.47 0.74 2.09 5,929.16 1.00
alpha[14] 0.01 0.61 −1.19 0.01 1.19 5,448.58 1.00
alpha[15] 1.71 0.77 0.33 1.66 3.40 5,404.62 1.00
alpha[16] 1.73 0.77 0.32 1.69 3.30 4,837.64 1.00
alpha[17] 2.56 0.70 1.34 2.51 4.08 4,757.44 1.00
alpha[18] 2.15 0.60 1.09 2.11 3.45 4,079.94 1.00
alpha[19] 1.81 0.55 0.81 1.78 2.98 6,034.53 1.00
alpha[20] 3.09 0.80 1.72 3.05 4.78 4,849.12 1.00
alpha[21] 2.16 0.61 1.08 2.12 3.49 4,492.36 1.00
alpha[22] 2.13 0.60 1.05 2.11 3.39 5,183.80 1.00
alpha[23] 2.15 0.62 1.03 2.11 3.45 5,733.02 1.00
alpha[24] 1.55 0.50 0.65 1.52 2.58 4,548.88 1.00
alpha[25] −1.10 0.44 −2.01 −1.09 −0.28 4,604.75 1.00
alpha[26] 0.09 0.38 −0.66 0.09 0.82 5,601.76 1.00
alpha[27] −1.54 0.51 −2.59 −1.51 −0.60 5,589.68 1.00
alpha[28] −0.56 0.40 −1.35 −0.55 0.19 5,049.25 1.00
alpha[29] 0.08 0.39 −0.66 0.08 0.86 5,158.46 1.00
alpha[30] 1.32 0.47 0.46 1.29 2.27 4,620.21 1.00
alpha[31] −0.72 0.41 −1.57 −0.71 0.09 5,258.36 1.00
alpha[32] −0.39 0.40 −1.24 −0.39 0.37 5,528.52 1.00
alpha[33] 2.83 0.67 1.68 2.77 4.24 4,391.22 1.00
alpha[34] 2.47 0.58 1.45 2.43 3.72 5,907.56 1.00
alpha[35] 2.48 0.61 1.40 2.44 3.81 3,674.38 1.00
alpha[36] 1.90 0.47 1.05 1.88 2.91 4,599.25 1.00
alpha[37] 1.92 0.49 1.04 1.90 2.91 4,875.96 1.00
alpha[38] 3.35 0.77 2.00 3.29 4.99 4,849.03 1.00
alpha[39] 2.46 0.58 1.44 2.43 3.71 5,110.38 1.00
alpha[40] 2.15 0.51 1.25 2.13 3.24 4,697.76 1.00
alpha[41] −1.91 0.50 −2.97 −1.89 −0.97 5,483.62 1.00
alpha[42] −0.63 0.36 −1.32 −0.62 0.06 5,997.09 1.00
alpha[43] −0.51 0.35 −1.21 −0.51 0.17 6,252.09 1.00
alpha[44] −0.40 0.34 −1.06 −0.39 0.24 4,936.93 1.00
alpha[45] 0.51 0.35 −0.16 0.50 1.21 4,893.38 1.00
alpha[46] −0.63 0.35 −1.36 −0.63 0.05 5,426.84 1.00
alpha[47] 1.91 0.49 1.00 1.88 2.92 5,090.97 1.00
alpha[48] −0.06 0.34 −0.72 −0.06 0.60 5,243.91 1.00
lp__ −524.64 4.91 −535.11 −524.27 −515.99 1,774.83 1.00

Podemos cotejar nuestras predicciones de este modelo No-pooling con los datos originales:

Predicciones del modelo 2 (no pooling)

Al evaluar las estimaciones del modelo No- pooling (un intercepto \(\alpha_i\) independiente por tanque) obtenemos:

  1. Gran variabilidad entre tanques
    • Los \(\alpha_i\) (en escala log-odds) oscilan aproximadamente entre \(-1.9\) y \(+3.4\).
    • Transformados a probabilidad, algunos tanques se estiman con supervivencias cercanas al 10–20% y otros al 90–95%.
  2. Incertidumbre muy desigual
    • Tanques con pocas observaciones (densidad pequeña) presentan desviaciones estándar de \(\alpha_i\) de 0.6–0.8 y rangos de credibilidad muy amplios (p. ej. \(\alpha_3: [–0.41, +2.06]\).
    • Tanques con más renacuajos reducen su incertidumbre a 0.3–0.5 en la desviación estándar de \(\alpha_i\).
  3. Sobreajuste
    • Las predicciones del modelo (puntos rojos) siguen casi exactamente los datos observados (puntos grises), incluso en valores extremos.
    • No existe “arrastre” hacia un promedio general: cada tanque se ajusta únicamente con su propia información.
  4. Problemas en tanques pequeños
    • Con muestras muy pequeñas, pocas muertes o supervivencias cambian drásticamente la estimación de \(\alpha_i\).
    • El ancho de los intervalos de credibilidad hace poco útiles esas predicciones para la toma de decisiones.

Por lo anterior, el modelo No-Pooling captura fielmente cada dato empírico, pero padece de sobreajuste y de alta incertidumbre en tanques con pocas observaciones. Para obtener estimaciones más estables y evitar extremos sin fundamento, es recomendable utilizar un modelo jerárquico (Partial Pooling) que comparta información entre tanques.

3. Modelo Partial Pooling

En este tercer modelo parcialmente agrupado (partial pooling), asignamos un intercepto distinto a cada tanque, pero también estimamos la variabilidad entre ellos. Esto nos permite captar la heterogeneidad entre tanques y, al mismo tiempo, compartir información entre ellos. Para esto, se agregan dos parámetros: \(\mu\) y \(\sigma\), los cuáles llamaremos hiperparámetros desde este punto.El punto es que, en el modelo anterior, todas las \(\alpha_i\) se distribuían Normal con una media y desviación estándar establecida o fija. Con esta modificación, los \(\alpha_i\) comparten una misma distribución, lo que le permite transmitir información entre tanques al modelo.

\[ S_i \sim \textrm{Binomial}(D_i,p_i) \]

\[ \textrm{logit}(p_i) = \alpha_{T[i]} \]

\[ \alpha_j = \textrm{Normal}(\mu, \sigma) \] \[ \mu \sim \textrm{Normal}(0, 1.5) \] \[ \sigma \sim \textrm{Exponential}(1) \]

Lo cual se traduce en el siguiente código de Stan:

data {
  int<lower=0> T;         // Num de tanques
  int<lower=0> S[T];      // Num de renacuajos que sobrevivieron
  int<lower=0> D[T];      // Densidad inicial
}

parameters {
  real<lower=0> mu_alpha;        // Promedio del alpha
  real<lower=0> sigma_alpha;      // Desv est de los alphas
  vector[T] alpha_tank;           // Alpha de cada tanque
}

model {
  // Hyperpriors 
  mu_alpha ~ normal(0, 1.5);
  sigma_alpha ~ exponential(1);
  // Priors
  alpha_tank ~ normal(mu_alpha, sigma_alpha);

  for (t in 1:T) {
    S[t] ~ binomial(D[t], inv_logit(alpha_tank[t]));
  }
}

generated quantities {
  int S_rep[T]; 

  for (t in 1:T) {
    S_rep[t] = binomial_rng(D[t], inv_logit(alpha_tank[t]));
  }
}
parameter mean sd 2.5% 50% 97.5% n_eff Rhat
alpha_tank[1] 2.12 0.89 0.55 2.03 4.04 5,139.11 1.00
alpha_tank[2] 3.08 1.10 1.19 2.97 5.49 3,122.16 1.00
alpha_tank[3] 1.00 0.68 −0.24 0.96 2.43 5,151.81 1.00
alpha_tank[4] 3.10 1.12 1.18 3.03 5.45 4,620.60 1.00
alpha_tank[5] 2.15 0.90 0.60 2.08 4.18 4,476.80 1.00
alpha_tank[6] 2.14 0.87 0.60 2.08 4.08 4,949.18 1.00
alpha_tank[7] 3.09 1.11 1.21 2.97 5.60 4,031.84 1.00
alpha_tank[8] 2.14 0.88 0.57 2.08 4.07 4,796.29 1.00
alpha_tank[9] −0.18 0.61 −1.40 −0.18 1.01 4,733.70 1.00
alpha_tank[10] 2.14 0.88 0.60 2.07 4.10 4,462.08 1.00
alpha_tank[11] 1.02 0.68 −0.21 0.98 2.44 5,086.39 1.00
alpha_tank[12] 0.59 0.63 −0.63 0.57 1.84 4,773.67 1.00
alpha_tank[13] 1.00 0.68 −0.27 0.98 2.39 5,379.87 1.00
alpha_tank[14] 0.20 0.60 −0.96 0.19 1.36 5,045.57 1.00
alpha_tank[15] 2.14 0.86 0.64 2.08 4.00 4,734.10 1.00
alpha_tank[16] 2.12 0.86 0.62 2.07 3.97 4,767.58 1.00
alpha_tank[17] 2.90 0.78 1.55 2.84 4.54 4,366.18 1.00
alpha_tank[18] 2.39 0.65 1.26 2.35 3.77 4,345.57 1.00
alpha_tank[19] 2.02 0.58 0.98 1.98 3.28 4,674.03 1.00
alpha_tank[20] 3.66 0.98 2.02 3.56 5.87 3,902.58 1.00
alpha_tank[21] 2.41 0.67 1.25 2.36 3.90 4,419.92 1.00
alpha_tank[22] 2.39 0.67 1.21 2.34 3.79 4,914.81 1.00
alpha_tank[23] 2.41 0.67 1.24 2.36 3.84 4,135.74 1.00
alpha_tank[24] 1.70 0.51 0.79 1.66 2.78 5,050.46 1.00
alpha_tank[25] −1.01 0.45 −1.93 −1.00 −0.17 5,079.17 1.00
alpha_tank[26] 0.16 0.39 −0.59 0.16 0.94 4,803.57 1.00
alpha_tank[27] −1.43 0.49 −2.45 −1.39 −0.52 4,242.92 1.00
alpha_tank[28] −0.46 0.41 −1.29 −0.46 0.32 4,962.52 1.00
alpha_tank[29] 0.17 0.41 −0.62 0.15 1.00 4,790.05 1.00
alpha_tank[30] 1.45 0.50 0.53 1.43 2.50 5,307.24 1.00
alpha_tank[31] −0.63 0.41 −1.47 −0.63 0.15 5,032.40 1.00
alpha_tank[32] −0.31 0.41 −1.12 −0.30 0.48 5,829.77 1.00
alpha_tank[33] 3.17 0.76 1.89 3.10 4.93 4,176.24 1.00
alpha_tank[34] 2.72 0.65 1.60 2.67 4.15 4,461.70 1.00
alpha_tank[35] 2.71 0.64 1.60 2.67 4.09 4,453.75 1.00
alpha_tank[36] 2.06 0.49 1.16 2.04 3.08 5,549.26 1.00
alpha_tank[37] 2.07 0.51 1.16 2.03 3.10 4,934.79 1.00
alpha_tank[38] 3.94 0.99 2.30 3.82 6.12 3,588.83 1.00
alpha_tank[39] 2.70 0.64 1.59 2.66 4.09 4,422.90 1.00
alpha_tank[40] 2.34 0.55 1.37 2.31 3.51 4,451.72 1.00
alpha_tank[41] −1.80 0.47 −2.77 −1.78 −0.96 4,372.33 1.00
alpha_tank[42] −0.57 0.36 −1.28 −0.57 0.11 5,267.18 1.00
alpha_tank[43] −0.45 0.35 −1.14 −0.45 0.21 5,703.00 1.00
alpha_tank[44] −0.34 0.33 −1.01 −0.34 0.30 3,975.99 1.00
alpha_tank[45] 0.58 0.35 −0.10 0.58 1.28 6,011.52 1.00
alpha_tank[46] −0.58 0.34 −1.26 −0.57 0.07 5,916.07 1.00
alpha_tank[47] 2.06 0.50 1.15 2.05 3.06 5,208.33 1.00
alpha_tank[48] 0.00 0.34 −0.66 0.00 0.67 4,967.75 1.00
lp__ −532.80 5.31 −544.04 −532.48 −523.18 1,616.43 1.00

Y para los hiperparámetros \(\mu\) y \(\sigma\):

parameter mean sd 2.5% 50% 97.5% n_eff Rhat
mu_alpha 1.35 0.26 0.85 1.34 1.85 2,947.94 1.00
sigma_alpha 1.62 0.21 1.26 1.60 2.10 2,549.44 1.00

Podemos nuevamente, ver cómo se comportan nuestras predicciones contra los valores observados:

Predicciones del modelo 3 (partial pooling)

Esta estructura jerárquica permite al modelo “aprender entre tanques”: los tanques con poca información se benefician del conocimiento global, mientras que los tanques con muchos datos retienen sus características particulares. Esto se conoce como regularización adaptativa o “shrinkage”.

Las estimaciones de los hiperparámetros:

  • La media global \(\mu_\alpha\) ≈ 1.35, lo que se traduce en una tasa de supervivencia promedio de aproximadamente 79% (vía función logística), es decir por arriba del promedio.

  • La desviación estándar \(\sigma_\alpha\) ≈ 1.62, lo cual indica una heterogeneidad considerable entre tanques.

Los interceptos de tanques con pocos renacuajos (menos información) se encogen más hacia la media global, generando estimaciones más conservadoras. En contraste, los tanques con mayor número de observaciones muestran menos shrinkage, permitiendo conservar la variación real.

El modelo Parcial Pooling, corrige el subajuste del modelo Total Pooling y el sobreajuste del modelo No-Pooling, logrando un balance óptimo.Como podemos ver en la última gráfica, las predicciones siguen de cerca los datos observados, pero sin reaccionar de forma extrema a variaciones pequeñas o poco informativas. También podemos observar en la tabla de resultados, que los valores de \(\hat{R}\) ≈ 1.00 y los tamaños efectivos (n_eff) son altos, lo cual indica que las cadenas MCMC convergieron adecuadamente y los resultados son confiables.

Con esto, si podemos concluir que el modelo jerárquico permite identificar patrones globales sin ignorar las diferencias individuales, siendo especialmente útil cuando algunos grupos (tanques) tienen poca información, con esto confirmas que existe variabilidad real en las tasas de supervivencia entre tanques.

Relación entre log-odds y probabilidad

Para comprender mejor el efecto del parámetro α en los modelos logísticos, se muestra a continuación la transformación inv_logit(α):

Esto ilustra cómo pequeñas diferencias en log-odds (α), se traducen en cambios más o menos pronunciados en la probabilidad, dependiendo de la región de la curva sigmoide.

Comparativos

Podemos graficar las estimaciones del modelo No-Pooling y el Parcial Pooling para ver el comportamiento con mayor detenimiento.

Incluyendo intervalos de credibilidad para ver su comportamiento:

Ahora que tenemos los tres modelos ajustados, es momento de compararlos. Para ello, utilizaremos:

  • Diagnósticos MCMC: Para confirmar la buena convergencia de los modelos.

  • Visualización de predicciones: Para contrastar la calidad de las predicciones por tanque.

Diagnósticos MCMC

Una buena práctica es comparar las trazas y distribuciones posteriores de parámetros clave para verificar convergencia y comportamiento estable.

Estos gráficos nos permiten confirmar que las cadenas exploran bien el espacio posterior, sin signos de problemas como falta de mezcla o divergencias.

Predicciones por tanque

Unificamos en una sola gráfica las predicciones de los tres modelos para compararlas directamente.

Este gráfico muestra cómo cada modelo aproxima los datos observados. Veremos que el Modelo Parcial Pooling logra un buen balance, sin subestimar ni sobreajustar.

Modelo con información completa

En los tres modelos anteriores, no incorporamos todas las variables, y como lo que se desea es estudiar el fenómeno, incluiyendo todas las variables, es decir, agregando la variable depredador P y la variable tamaño G, dado nuestro DAG, puede mejorar la estimación de los efectos por variable. Podemos extender entonces la idea del Modelo Jerárquico e incorporar estas variables. De tal modo que:

\[ S_i \sim \textrm{Binomial}(D_i, p_i) \]

\[ \textrm{logit}(p_i) = \alpha_{\textrm{Tanque}[i]} + \beta_p *\textrm{pred} + \beta_s *\textrm{size}\]

\[ \alpha_j \sim \textrm{Normal}(\mu, \sigma) \]

\[ \mu \sim \textrm{Normal}(0, 1.5) \]

\[ \sigma \sim \textrm{Exponential}(1) \]

\[ \beta_p \sim \textrm{Normal}(-0.5,1) \]

\[ \beta_s \sim \textrm{Normal}(0,1) \]

Por lo que el codigo en stan:

data {
  int<lower=0> N;               // Número total de observaciones
  int<lower=0> T;               // Número de tanques
  int<lower=0> surv[N];         // renacuajos que sobrevivieron
  int<lower=0> density[N];      // Densidad inicial de renacuajos
  int<lower=1, upper=T> tank[N];// Índice del número de tanque
  vector[N] pred;               // Presencia de depredadores (0 o 1)
  vector[N] size;               // Tamaño de los renacuajos
}

parameters {
  real<lower=0> mu_alpha;      // Promedio de alpha
  real<lower=0> sigma_alpha;    // Desviación estándar de los alphas
  vector[T] alpha_tank;         // Alpha para cada tanque
  real bp;                      // Coef para los depredadores
  real bs;                      // Coef del tamaño de los renacuajos
}

model {
  // Hyperpriors
  mu_alpha ~ normal(0, 1.5);
  sigma_alpha ~ exponential(1);

  // Priors
  alpha_tank ~ normal(mu_alpha, sigma_alpha);
  bp ~ normal(-0.5, 1);
  bs ~ normal(0, 1);  // Prior normal para el coeficiente de tamaño

  // Modelo para cada observación
  for (n in 1:N) {
    real p = inv_logit(alpha_tank[tank[n]] + bp * pred[n] + bs * size[n]);
    surv[n] ~ binomial(density[n], p);
  }
}

generated quantities {
  int surv_rep[N]; // Predicciones replicadas para la supervivencia

  for (n in 1:N) {
    surv_rep[n] = binomial_rng(density[n], inv_logit(alpha_tank[tank[n]] + bp * pred[n] + bs * size[n]));
  }
}
parameter mean sd 2.5% 50% 97.5% n_eff Rhat
alpha_tank[1] 2.76 0.70 1.41 2.72 4.16 2,434.28 1.00
alpha_tank[2] 3.17 0.72 1.83 3.14 4.68 2,503.35 1.00
alpha_tank[3] 2.03 0.64 0.76 2.04 3.28 2,016.90 1.00
alpha_tank[4] 3.17 0.73 1.80 3.16 4.65 2,635.81 1.00
alpha_tank[5] 2.63 0.69 1.30 2.62 4.06 2,491.62 1.00
alpha_tank[6] 2.63 0.68 1.32 2.61 4.04 3,909.75 1.00
alpha_tank[7] 3.09 0.74 1.75 3.05 4.62 3,750.74 1.00
alpha_tank[8] 2.62 0.68 1.34 2.59 3.98 3,671.17 1.00
alpha_tank[9] 2.55 0.62 1.30 2.57 3.72 1,106.12 1.00
alpha_tank[10] 3.83 0.63 2.62 3.81 5.10 1,159.63 1.00
alpha_tank[11] 3.30 0.63 2.08 3.28 4.55 1,184.57 1.00
alpha_tank[12] 3.04 0.60 1.85 3.04 4.21 1,119.06 1.00
alpha_tank[13] 3.08 0.59 1.94 3.06 4.28 1,566.93 1.00
alpha_tank[14] 2.57 0.56 1.47 2.56 3.68 1,851.18 1.00
alpha_tank[15] 3.63 0.62 2.45 3.62 4.91 1,715.61 1.00
alpha_tank[16] 3.61 0.61 2.45 3.59 4.90 2,085.91 1.00
alpha_tank[17] 3.16 0.62 1.98 3.13 4.41 2,343.04 1.00
alpha_tank[18] 2.85 0.59 1.76 2.84 4.09 2,274.04 1.00
alpha_tank[19] 2.60 0.56 1.51 2.59 3.73 2,052.29 1.00
alpha_tank[20] 3.51 0.69 2.27 3.47 4.95 2,480.86 1.00
alpha_tank[21] 2.65 0.55 1.61 2.63 3.77 4,021.02 1.00
alpha_tank[22] 2.65 0.57 1.57 2.63 3.82 4,108.91 1.00
alpha_tank[23] 2.66 0.55 1.64 2.63 3.81 4,017.67 1.00
alpha_tank[24] 2.11 0.50 1.18 2.08 3.14 3,521.85 1.00
alpha_tank[25] 1.94 0.55 0.84 1.96 2.96 777.49 1.00
alpha_tank[26] 2.89 0.52 1.84 2.90 3.89 682.28 1.01
alpha_tank[27] 1.63 0.58 0.42 1.66 2.71 824.78 1.00
alpha_tank[28] 2.37 0.52 1.34 2.36 3.36 695.80 1.01
alpha_tank[29] 2.58 0.44 1.70 2.60 3.46 1,064.99 1.00
alpha_tank[30] 3.54 0.49 2.62 3.54 4.50 1,174.95 1.00
alpha_tank[31] 1.93 0.48 0.95 1.94 2.83 1,013.02 1.00
alpha_tank[32] 2.19 0.46 1.25 2.20 3.06 996.11 1.00
alpha_tank[33] 3.34 0.60 2.21 3.32 4.59 2,464.18 1.00
alpha_tank[34] 3.06 0.56 2.00 3.06 4.23 2,098.78 1.00
alpha_tank[35] 3.06 0.58 2.00 3.04 4.26 2,422.46 1.00
alpha_tank[36] 2.60 0.52 1.63 2.59 3.62 1,795.38 1.00
alpha_tank[37] 2.34 0.49 1.43 2.32 3.38 3,364.41 1.00
alpha_tank[38] 3.49 0.64 2.37 3.45 4.86 3,015.16 1.00
alpha_tank[39] 2.83 0.53 1.84 2.82 3.95 4,604.74 1.00
alpha_tank[40] 2.58 0.52 1.64 2.56 3.65 3,765.14 1.00
alpha_tank[41] 1.31 0.57 0.11 1.33 2.36 738.26 1.00
alpha_tank[42] 2.28 0.51 1.23 2.30 3.23 657.09 1.01
alpha_tank[43] 2.39 0.51 1.33 2.40 3.36 669.88 1.01
alpha_tank[44] 2.49 0.50 1.42 2.50 3.46 666.02 1.01
alpha_tank[45] 2.94 0.42 2.10 2.93 3.76 953.41 1.00
alpha_tank[46] 1.95 0.44 1.04 1.97 2.78 1,017.15 1.00
alpha_tank[47] 4.02 0.49 3.09 4.00 5.06 1,287.33 1.00
alpha_tank[48] 2.44 0.42 1.60 2.44 3.26 944.36 1.00
lp__ −501.01 6.88 −515.32 −500.76 −488.35 924.49 1.00

Y para los hiperparámetros \(\mu\), \(\sigma\) y las \(\beta\)´s :

parameter mean sd 2.5% 50% 97.5% n_eff Rhat
mu_alpha 2.74 0.28 2.19 2.74 3.27 484.18 1.01
sigma_alpha 0.78 0.14 0.53 0.78 1.09 1,214.49 1.00
bp −2.44 0.30 −3.01 −2.45 −1.84 567.64 1.01
bs −0.41 0.29 −0.95 −0.41 0.18 887.21 1.00

Como en los modelos anteriores, podemos ver cómo se comportan nuestras predicciones de este modelo completo contra los valores observados:

Si comparamos este modelo completo contra el Parcial Pooling que era hasta este momento el mejor modelo, podemos visualizar esta comparación:

Al incluir las variables de presencia de depredadores y tamaño, el Modelo Completo redistribuye parte de la variabilidad extrema hacia efectos sistemáticos: las predicciones de supervivencia caen en los tanques con depredadores y suben en los tanques sin ellos, tal como cabe esperar desde el punto de vista biológico.

Esa mayor separación aparente en la zona de altas tasas de supervivencia, no indica un error, sino que el modelo está capturando la variación mediante las covariables y suavizando (regularizando) los valores extremos.

Conclusiones

En este proyecto se implementaron modelos jerárquicos bayesianos para analizar datos experimentales sobre la supervivencia de renacuajos, recolectados por Vonesh & Bolker en un entorno controlado con variaciones en densidad poblacional y presencia de depredadores. A través de un enfoque comparativo, se evaluó tres estrategias de modelación —total pooling, no pooling y partial pooling— y se exploró sus implicaciones estadísticas y biológicas.

Los principales hallazgos son:

Los modelos Total Pooling subestiman la variabilidad

El enfoque de total pooling ignora la estructura jerárquica de los datos, asumiendo una única probabilidad de supervivencia para todos los tanques. Esto conduce a un subajuste (underfitting), perdiendo de vista diferencias importantes debidas a factores no observados, como la variación ambiental o genética entre tanques.

Los modelos No-pooling sobreajustan

El modelo sin agrupamiento (no pooling), asigna parámetros independientes a cada tanque, sin compartir información entre ellos. Esto resulta en sobreajuste (overfitting), especialmente en tanques con pocos datos, donde incluso pequeñas fluctuaciones generan estimaciones inestables.

El modelo Jerárquico o Parcial Pooling logra un equilibrio

El modelo parcialmente agrupado (partial pooling), combina la flexibilidad de asignar parámetros específicos a cada grupo con la estabilidad de una distribución común que los regula. Esta estructura:

  • Captura la heterogeneidad real entre tanques.

  • Reduce la incertidumbre en grupos pequeños mediante shrinkage adaptativo.

  • Mejora la capacidad predictiva y generalización del modelo.

A nivel biológico, se encontró que:

  • La supervivencia promedio fue mayor en renacuajos pequeños, un resultado que, aunque contraintuitivo, coincide con las conclusiones del estudio original.
  • La densidad inicial afecta negativamente la supervivencia, probablemente por competencia intraespecífica o mayor exposición a depredadores.
  • La presencia de depredadores incrementa la variabilidad de las tasas de supervivencia, posiblemente por efectos de habituación o respuestas defensivas plásticas.
  • Las distribuciones iniciales de supervivencia pueden actuar como regularizadores informativos, permitiendo una inferencia más robusta.

Notas

  1. Vonesh, J. R., & Bolker, B. M. (2005). Statistical tools for analyzing larval amphibian survival data. Ecology, 86(1), 172-182.↩︎

  2. Se refiere a la capacidad de los embriones de la rana para ajustar su desarrollo y eclosión en respuesta a cambios ambientales, como la presencia de depredadores o el secado de su hábitat↩︎